%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 6Li D2 line, yet to consider are Doppler shifts.
%
% M: Complete Bloch matrix for the 6Li-D2-System.
% r: [ r11,r22,..,rii,Im(r12),Re(r12),Im(r13),Re(r13),..,Im(rhi),Re(rhi) ]
% b: Inhomogeneities only for the trace-reduced density matrix
%
%			mj=-5/2..+5/2	--- |d> --- |e> --- |f> --- |g> --- |h> --- |i>
%			mj=-3/2..+3/2			--- |9> --- |a> --- |b>	--- |c>
%	P3/2	mj=-1/2..+1/2					--- |7> --- |8>
%
%									 
%			mj=-3/2..+3/2			--- |3> --- |4> --- |5> --- |6>
%	S1/2	mj=-1/2..+1/2					--- |1> --- |2>
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function M = Lithium6LevelsTime
global Lg Le Jg Je S I K hf1 hf2 hf3

Lg =  0;		Le =  1;														%orbital angular momenta
Jg = .5;		Je =1.5;														%total angular momenta
S  = .5;		I  =  1;		K =  1;											%electron spin  &  nuclear spin  &  dipolar order
hf1=-2.8e6;		hf2=  0;		hf3=1.7e6;										%hyperfine structure

Opb=    3e6;					Orp=   .2e6;									%Rabi frequencies
Dpb=	  0;					Drp=      0;									%detunings for probe and repumper
Epb=[1 i 0];			        Erp=[1 i 0];									%polarisation vectors
Lpb=    1e6;					Lrp=    1e6;									%laser linewidths
B  =[0 0 1]'*1e-4;																%B-field
M  =Lithium6LevelsTimeMatrix(Opb,Orp,Dpb,Drp,Epb,Erp,Lpb,Lrp,B);
r0 =zeros(324,1);				r0(3:6)=1/4;									%initial populations
t  =[0:.05:1]*20e-6;

% fig=figure;
% set(1,'DoubleBuffer','on');
% mov=avifile('C:\Dokumente und Einstellungen\Philippe\Eigene Dateien\Arbeit\example3.avi','Quality',100,'FPS',1);
for k=1:length(t)
	r=expm(M*t(k))*r0;
	R(1:18,k)=r(1:18);
	PlotPops(R,t(1:k));
	h=PlotDots(r,t(k));
% 	set(h,'EraseMode','xor');
% 	F  =getframe(gcf);
% 	mov=addframe(mov,F);
end
% mov = close(mov);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting the populations

function h = PlotDots(r,t)

X=[1 1];						L=[-.3 .3];
figure(2);						clf;
for kk=1:18
    [Fn,Mn]=Q(kk);              energ=Fn/5*(1-.5*(kk>6))+.8*(kk>6)-.2;
    txt1=['h=plot(',num2str(Mn),'+L,',num2str(energ),'*X,''k'',',num2str(Mn),',',num2str(energ),',''or'','];
    txt2=['''LineWidth'',1.5,''Markersize'',r(',num2str(kk),')*50+1,''MarkerFaceColor'',''r'');'];
    eval([txt1,txt2]);          hold on;
end
set(gca,'xLim',[-3 3],'yLim',[-.1 1.1],'FontSize',16,'LineWidth',1);
set(gcf,'Color','w');
xLabel('\fontsize{20}{\itm_F}');
ylabel('\fontsize{20}{\itE}');
title(['\fontsize{16}{\itt} = ',num2str(t*1e6),' \mus']);
axis off;
pause(.01);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Plotting the populations

function PlotPops(R,t)

figure(1);						clf;
subplot(331),plot(t*1e6,R(7,:),'r',t*1e6,R(8,:),'b');
set(gca,'yLim',[0 1],'FontSize',16,'LineWidth',1);
ylabel('\fontsize{16}|{\itF''}=1/2\rangle');
subplot(332),plot(t*1e6,R(9,:),'c',t*1e6,R(10,:),'r',t*1e6,R(11,:),'b',t*1e6,R(12,:),'m');
set(gca,'yLim',[0 1],'FontSize',16,'LineWidth',1);
ylabel('\fontsize{16}|{\itF''}=3/2\rangle');
subplot(333),plot(t*1e6,R(13,:),'g',t*1e6,R(14,:),'c',t*1e6,R(15,:),'r',t*1e6,R(16,:),'b',t*1e6,R(17,:),'m',t*1e6,R(18,:),'k');
set(gca,'yLim',[0 1],'FontSize',16,'LineWidth',1);
ylabel('\fontsize{16}|{\itF''}=5/2\rangle');
subplot(334),axis off;
text(.1,1.,'\fontsize{12}green   {\itm}=-5/2');
text(.1,.8,'\fontsize{12}cyan    {\itm}=-3/2');
text(.1,.6,'\fontsize{12}red     {\itm}=-1/2');
text(.1,.4,'\fontsize{12}blue    {\itm}=+1/2');
text(.1,.2,'\fontsize{12}magenta {\itm}=+3/2')
text(.1,.0,'\fontsize{12}black   {\itm}=+5/2');
subplot(337),plot(t*1e6,R(1,:),'r',t*1e6,R(2,:),'b');
set(gca,'yLim',[0 1],'FontSize',16,'LineWidth',1);
xlabel('\fontsize{16}{\itt}   (\mus)');
ylabel('\fontsize{16}|{\itF}=1/2\rangle');
subplot(339),plot(t*1e6,R(3,:),'c',t*1e6,R(4,:),'r',t*1e6,R(5,:),'b',t*1e6,R(6,:),'m');
set(gca,'yLim',[0 1],'FontSize',16,'LineWidth',1);
xlabel('\fontsize{16}{\itt}   (\mus)');
ylabel('\fontsize{16}|{\itF}=3/2\rangle');
pause(.01);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Submatrices of DETUNING:
%    	Im r_mn		Re r_mn
%  (	  0			+-Dmn	)	Im r_mn
%  (	+-Dmn		  0		)	Re r_mn
%  The detunings are (0 -Dmn;+Dmn 0) if the energy increases from |m> to |n> and inverted else.
%
% HFS- and ZEEMAN-splitting:
%  to the (pos.) detuning the energy shift is added inwards and subtracted outwards
%
% Submatrices COHERENCE - COHERENCE:
%    	Im r_ab		Re r_ab
%  (	  0			+-Wmn	)	Im r_xy
%  (	+-Wmn		  0		)	Re r_xy
%  Above the diagonal, the Rabi frequencies are (0 -Wby;+Wby 0) if a==x or (0 -Wby;-Wby 0) if a==y or (0 -Wby;+Wby 0) if b==y 
%   or (0 +Wby;+Wby 0) if b==x
%
% Submatrices POPULATION - COHERENCE:
%  Connects ...
%
% For all of the above components, the Bloch matrix is antisymmetric. It is not for the following components.
%
% Submatrices DECAY of COHERENCES:
%  Every coherence between a pair of states decays with half the sum of all decay channels of both states.
%  For those coherences which are directly laser-driven, the laser linewidth is added.
%
% Submatrices DECAY of POPULATION:
%  In the diagonal are the decay width of the channels.
%  The decays must all be compensated by gain in other channels.

function M = Lithium6LevelsTimeMatrix(Opb,Orp,Dpb,Drp,Epb,Erp,Lpb,Lrp,B)
G  =6e6;																		%natural decay

M=zeros(324,324);
for ee=7:18																		%decay of excited states
	M(ee,ee)=-G;
end
for gg=1:6																		%filling of ground states
	for ee=7:18	
		M(gg,ee)=G * Transition(gg,ee);
	end
end
for gg=1:18																		%decay of coherences
	for ee=7:18
		M(C(gg,ee)  ,C(gg,ee)  )=-G/2 - G/2*(gg>=7) + Lpb*(gg>2&gg<7) + Lrp*(gg<=2);
		M(C(gg,ee)+1,C(gg,ee)+1)=-G/2 - G/2*(gg>=7) + Lpb*(gg>2&gg<7) + Lrp*(gg<=2);
	end
end

N=zeros(324,324);
for gg=1:6																		%population - coherences
	for ee=7:18
		N(gg,C(gg,ee))= Rabi(gg,ee,Opb,Orp,Epb,Erp,B) * Transition(gg,ee);
		N(ee,C(gg,ee))=-Rabi(gg,ee,Opb,Orp,Epb,Erp,B) * Transition(gg,ee);
	end
end
for gg=1:6																		%coherences - coherences
	for ee=7:18
		for ff=1:18
			if ff~=gg & ff~=ee
				sgA=(gg==min(ff,ee) | ee==max(ff,ee))*2-1;
				sgB=(gg==max(ff,ee) | ee==min(ff,ee))*2-1;
				N(C(gg,ff),C(ff,ee)+1)= sgA * Rabi(gg,ee,Opb,Orp,Epb,Erp,B) * Transition(gg,ee) /2;
				N(C(gg,ff)+1,C(ff,ee))= sgB * Rabi(gg,ee,Opb,Orp,Epb,Erp,B) * Transition(gg,ee) /2;
			end
		end
	end
end
for gg=1:6																		%detunings with upper hfs splitting & Zeeman shift
	for ee=7:18
		N(C(gg,ee),C(gg,ee)+1)=-Dpb*(gg<=2) - Drp*(gg>=3) + Zeeman(gg,ee,B);
	end
end
for gg=1:2																		%combined detunings
	for ff=3:6
		if gg~=ff
			N(C(gg,ff),C(gg,ff)+1)=-Dpb + Drp + Zeeman(gg,ff,B) - Zeeman(gg,ff,B);
		end
	end
end
M=M+N-N';


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Index subroutine

function ind = C(n1,n2)

m1=min(n1,n2);				m2=max(n1,n2);
vec=[];
for kk=1:18 vec=[vec,kk+[kk+1:18]/100]; end
dum=17+2*[1:length(vec)];
ind=dum(m1+m2/100==vec);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Quantum number subroutine

function [F,M] = Q(nn)
global Jg Je I

Fm=[];                      Mm=[];
for f=[abs(Jg-I):Jg+I abs(Je-I):Je+I]
    Fm=[Fm f*ones(1,2*f+1)];
    Mm=[Mm -f:f];
end
S =struct('Level',[1:18],'F',Fm,'M',Mm);
F =S.F(nn);					M =S.M(nn);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rabi frequency selection subroutine, old method:
%  Www = acos(Eee*B/norm(Eee)/norm(B))
%  Sss=Ooo*cos(Www)*(Me==Mg) + Ooo*sin(Www)*(abs(Me-Mg)==1);

function Sss = Rabi(gg,ee,Opb,Orp,Epb,Erp,B)

Ooo=(Orp*(gg<=2) + Opb*(gg>=3)) * (ee>=7);
Eee=(Erp*(gg<=2) + Epb*(gg>=3)) * (ee>=7);
[Fg,Mg]=Q(gg);	[Fe,Me]=Q(ee);
Sss=Ooo*Hanle(Eee/norm(Eee),B,Me-Mg);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Linear or circular polarization for pi or sigma transitions.

function Ppp = Hanle(Eee,B,dM)

g =[0 1 0]';
e3=B/norm(B);
e2=Kreuz(e3,g )/norm(Kreuz(e3,g ));
e1=Kreuz(e2,e3)/norm(Kreuz(e2,e3));
ep=(e1-i*e2)/sqrt(2);
em=(e1+i*e2)/sqrt(2);
Ppp=(Eee*ep)^2*(dM==1) + (Eee*e3)^2*(dM==0) + (Eee*em)^2*(dM==-1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Vector product

function c = Kreuz(a,b)

c=[a(2)*b(3)-a(3)*b(2)
   a(3)*b(1)-a(1)*b(3)
   a(1)*b(2)-a(2)*b(1)];


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Zeeman shift

function Z = Zeeman(gg,ee,B)
Constant;
global Lg Le Jg Je hf1 hf2 hf3

[Fg,Mg]=Q(gg);	[Fe,Me]=Q(ee);
Z = hf1*(ee<=8) + hf2*(ee>8&ee<13) + hf3*(ee>=13) + uB*norm(B)/(2*pi*hP)*(Gfactor(Lg,Jg,Fg)*Mg - Gfactor(Le,Je,Fe)*Me);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% g-factors

function gF = Gfactor(L,J,F)
Constant;
global S I

gJ= 1+(J*(J+1)+S*(S+1)-L*(L+1))/2/J/(J+1);
gF=gJ*(F*(F+1)+J*(J+1)-I*(I+1))/2/F/(F+1) + gJ*uK/uB*(F*(F+1)-J*(J+1)+I*(I+1))/2/F/(F+1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Rabi frequency between level |g> and |e>

function S = Transition(gg,ee)
global Jg Je I K

S =struct('Level',[1:18],'F',[.5*ones(1,2) 1.5*ones(1,4) .5*ones(1,2) 1.5*ones(1,4) 2.5*ones(1,6)],...
	'M',[-.5:.5 -1.5:1.5 -.5:.5 -1.5:1.5 -2.5:2.5]);
[Fg,Mg]=Q(gg);	[Fe,Me]=Q(ee);
S =Symbol3j(Fg,K,Fe,Mg,sign(Me-Mg),Me)^2*Symbol6j(Je,Jg,K,Fg,Fe,I)^2*(2*Fg+1)*(2*Je+1)*(2*K+1)/(2*I+1);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates the Clebsch-Gordon coefficient (j1 j2 j3)
%                      						{m1,m2,m3}

function C = Symbol3j(j1,j2,j3,m1,m2,m3)

if m3 ~= m1+m2, C=0; return; end
vmin = max([-j1+j2+m3,-j1+m1,0]);
vmax = min([j2+j3+m1,j3-j1+j2,j3+m3]);
C=sqrt((2*j3+1)*prod(1:j3+j1-j2)*prod(1:j3-j1+j2)*prod(1:j1+j2-j3)*prod(1:j3+m3)*prod(1:j3-m3)...
			/prod(1:j1+j2+j3+1)/prod(1:j1-m1)/prod(1:j1+m1)/prod(1:j2-m2)/prod(1:j2+m2));
S=0;
for v=vmin:vmax
   S=S+(-1)^(v+j2+m2)/prod(1:v)*prod(1:j2+j3+m1-v)*prod(1:j1-m1+v)/prod(1:j3-j1+j2-v)/prod(1:j3+m3-v)/prod(1:v+j1-j2-m3);
end
C=C*S;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates 6j-symbols {j1 j2 j3}
%                       {J1 J2 J3}

function sym = Symbol6j(j1,j2,j3,J1,J2,J3)

un=ceil(max([0 j1+j2+j3 j1+J2+J3 J1+j2+J3 J1+J2+j3]));
ob=floor(min([j1+j2+J1+J2 j2+j3+J2+J3 j1+j3+J1+J3]));
jsym=0;
for t=un:1:ob
	jsym=jsym+(-1)^t*prod(1:t+1)/prod(1:t-j1-j2-j3)/prod(1:t-j1-J2-J3)/prod(1:t-J1-j2-J3)/prod(1:t-J1-J2-j3)...
								/prod(1:j1+j2+J1+J2-t)/prod(1:j2+j3+J2+J3-t)/prod(1:j1+j3+J1+J3-t);
end
sym=sqrt(Racah(j1,j2,j3)*Racah(j1,J2,J3)*Racah(J1,j2,J3)*Racah(J1,J2,j3)) * jsym;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Calculates Racahs triangle coefficient.

function rac = Racah(a,b,c)

f1=prod(1:-a+b+c);
f2=prod(1: a-b+c);
f3=prod(1: a+b-c);
f4=prod(1: a+b+c+1);
rac=f1*f2*f3/f4;
